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t 

FOREWORD 


An  eight  node,  first  order  hexahedron  finite  element  has  been  modified  to 
simulate  an  ideal  compressible  fluid.  The  transient  response  of  two  submerged 
fluid-coupled  cylindrical  elastic  shells  impinged  upon  by  a  transversely  incident 
plane  pressure  pulse  is  computed  using  the  fluid  elements  to  model  the  fluid 
between  the  two  shells.  The  numerical  computations  are  accomplished  by  utlizing 
the  USA-STAGS  computer  code.  This  code  employs  the  finite  element  method  for  the 
structure  and  the  Doubly  Asymptotic  Approximation  for  the  fluid-structure  inter¬ 
action.  The  results  of  the  USA-STAGS  calculations  are  compared  to  an  analytical 
solution. 

This  work  was  completed  as  part  of  the  Double  Hull  Lethality  Program  and 
was  funded  by  the  Defense  Nuclear  Agency. 

The  author  acknowledges  the  help,  comments,  and  suggestions  provided  by 
Michael  Giltrud,  Frank  Brogan,  and  Dr.  Hanson  Huang. 
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INTRODUCTION 

The  class  of  problems  involving  fluid  structure  interaction  is  indeed  a  com¬ 
plex  one.  They  simultaneously  bring  together  the  subjects  of  both  fluid  mechanics 
and  solid  mechanics.  The  structure  must  first  be  modeled  adequately  to  determine 
the  elastic  response.  This  is  certainly  within  the  capability  of  many  structural 
analysis  computer  codes.  However,  for  a  structure  submerged  in  a  fluid,  the 
surface  loading  depends  on  the  state  of  motion  of  the  structure  in  the  fluid. 

Thus,  the  problem  of  determining  the  loading  is  one  in  which  the  state  of  motion 
of  the  fluid  and  structure  are  coupled. 

Analytical  investigations  have  resulted  in  the  adoption  of  methods  for 
uncoupling  the  structure  from  the  fluid.  Several  such  Drocedures,  known  as  sur¬ 
face  approximation  techniques,  have  been  compared  by  Geers1.  It  is  apparent  from 
his  analysis  that  the  Doubly  Asymptotic  Approximation  (DAA)1.2*3  is  the  most 
accurate  for  early  and  late  times  in  predicting  the  behavior  of  submerged  shells. 
The  DAA  also  affords  a  smooth  transition  between  early  and  late  time  response. 

The  STAGS  finite  element  code  (Structural  Analysis  of  General  Shells)4 
has  been  combined  with  a  DAA  code,  USA  (Underwater  Shock  Analysis)5  with  the 


Geers,  T.  L.,  "Transient  Response  Analysis  of  Submerged  Structures,"  in  Finite 
Element  Analysis  of  Transient  Non-Linear  Behavior,  AMD  Vol .  14,  ASME,  New  York, 
WIT. -  - 

2Geers,  T.  L.,  "Response  of  an  Elastic  Cylindrical  Shell  to  a  Transverse  Acoustic 
Shock  Wave  In  a  Liqht  Fluid  Medium,"  J.  Acoust.  Soc.  Am.,  Vol.  48,  No.  3,  Sep  1970, 
pp.  692-701. 

3Geers,  T.  1.,  "Excitation  of  an  Elastic  Cylindrical  Shell  by  a  Transient  Acoustic 
Wave,"  J.  Apol.  Mech.,  Vol.  36,  No.  3,  Sep  1969,  op.  459-469. 

4A1mroth,  B.  0.,  Brogan,  F.  A.,  "The  STAGS  Computer  Code,"  Report  No.  LF'’SC-DS58853, 
Structural  Mechanics  Laboratory,  Lockheed  Palo  Alto  Research  Lab,  Falo  Alto,  CA. 

5DeRuntz,  J.  A.,  Geers,  T.  L.,  Fell  1  pa,  C.  A.,  "The  Underwater  Shock  Analysis  (USA) 
Code,  A  Reference  Manual,"  LMSC-D624328,  Contract  No.  DNA  001-76-C-0285. 

28  Feb  1978. 
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resulting  code  called  USA-STAGS6.  This  code  has  been  used  to  predict  the  response 
of  a  linear  elastic  cylindrical  shell  immersed  in  a  fluid  through  which  a  plane 
acoustic  wave  is  propagating7.  The  results  of  the  USA-STAGS  calculations  have 
been  compared  to  an  exact  analysis  by  Huang8.  Good  agreement  has  been  obtained  for 
step  and  exponential  pressure  wave  loading.  In  addition,  the  USA-STAGS  predictions 
compare  favorably  with  the  late  time  asymptote. 

This  paper  investigates  the  effectiveness  of  using  a  fluid  finite  element  to 
model  the  entrained  fluid  between  two  submerged  concentric  cylindrical  elastic 
shells  impinged  upon  by  a  transversely  incident  pressure  pulse.  The  results  of 
the  USA-STAGS  calculations  are  compared  to  an  analytical  solution  by  Huang9. 

ANALYTICAL  FORMULATION 


The  USA-STAGS  code  is  a  combination  of  the  Underwater  Shock  Analysis  (USA) 
code  and  the  Structural  Analysis  of  General  Shells  code  (STAGS).  In  USA,  the 
fluid  is  assumed  to  be  an  infinite  acoustic  medium  whose  response  to  the  motions 
of  the  structure  is  described  by  the  DAA.  STAGS  is  a  general  purpose  finite 
element  code  Intended  for  analysis  of  a  shell  type  structure. 

At  this  time  it  is  useful  to  highlight  the  essentials  of  the  problem  formula¬ 
tion. 

STRUCTURAL  RESPONSE  EQUATION.  As  STAGS  is  based  upon  the  finite  element 
method,  the  discretized  differential  equation  of  motion  for  the  non-linear  struc¬ 
ture  is  expressed  as 


Msx  +  Csx  +  Ksx  *  f  (1) 

where  x  Is  the  structural  displacement  vector.  M$  and  Cg  are  the  structural  mass 
and  damping  matrices.  K$  Is  the  non-linear  stiffness  matrix  and  f  is  the  external 
force  vector.  Generally  M$,  C$  and  Ks  are  highly  banded  symmetric  matrices  of 
large  order.  In  particular,  STAGS  considers  Ms  to  be  diagonal  and  Cg  to  be  a 
linear  combination  of  Mg  and  !<s. 


sDeRuntz,  0.  A.,  Brogan,  F.  A.,  "Underwater  Shock  Analysis  of  Non-Linear  Structures, 
A  Reference  Manual  for  the  USA-STAGS  Code,"  Report  No.  LMSC-C633864,  Lockheed  Palo 
Alto  Research  Laboratory,  Palo  Alto,  CA,  Jan  1979. 

7Giltrud,  M.  E.,  and  Lucas,  D.  S.,  "A  Numerical  Comparison  with  an  Exact  Solution 
for  the  Transient  Response  of  a  Cylinder  Immersed  in  a  Fluid,"  in  the  49th  Shock 
and  Vibration  Bulletin,  Sep  1979. 

8Huang,  H.,  "An  Exact  Analysis  of  the  Transient  Interaction  of  Acoustic  Plane 
Waves  with  a  Cylindrical  Elastic  Shell,"  J.  Appl.  Mech.,  Vol .  37,  Mo.  4,  Dec  1970, 
pp.  1091- 11C6. 

9Huang,  H.,  "Transient  Response  of  Two  Fluid-Coupled  Cylindrical  Elastic  Shells  to 
an  Pressure  Pulse,"  ASME  Journal  of  Applied  Mechanics,  Vol.  46,  1979, 
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For  the  excitation  of  a  submerged  structure  by  a  transient  acoustic  wave,  f 
is  given  by 

f  -  -e  Af  (Bl  +  es)  +  rD  (2) 

where  p^  is  the  modal  incident  pressure  vector  (a  known)  and  p$  is  the  modal 
scattered  pressure  vector  (unknown).  The  dry  structure  dynamic  load  vector  is 
given  by  additionally,  Af  is  an  area  matrix  and  G  is  a  transformation  matrix. 

FLUID  RESPONSE  EQUATION.  USA  makes  use  of  the  DAA  to  describe  the  response 
of  the  scattered  pressure  at  the  fluid  structure  interface10*11.  The  DAA  exhibits 
both  excellent  high  frequency  accuracy  and  excellent  low  frequency  accuracy  as  well 
as  offering  a  smooth  transition  between  the  two  asymptotes. 

The  differential  equation  governing  the  fluid  response  is 

Mf  P$  +  pc  p$  =  pc  Mf  Us  (3) 

where  p$  is  the  scattered  pressure  vector;  ys  is  the  vector  of  the  scattered  wave 
particle  velocities;  p  and  c  are  the  fluid  density  and  sound  speed.  The  added 
mass  matrix,  Mf,  is  produced  by  a  boundary  -  element  treatment  of  Laplace's 
equation  for  the  irrotational  flow  generated  in  an  inviscid,  incompressible  fluid 
by  motions  of  the  structure's  wet  surface12. 

The  above  equation  (3)  is  subject  to  the  following  kinematic  compatibility 
equation 

GTx  =  Ur  +  Us  (4) 

where  the  superscript  T  represents  the  matrix  transposition.  The  compatibility 
equation  (4)  constrains  the  normal  fluid  particle  velocity  (Uj  +  U$)  to  the 
normal  structural  velocity  at  the  wet  interface.  The  transformation  matrix,  G, 
relates  the  structural  degrees  of  freedom  to  the  fluid  degrees  of  freedom  and  it 
follows  from  the  invariance  of  virtual  work  with  respect  to  coordinate  systems. 


10Geers,  T.  L.,  "Residual  Potential  and  Approximate  Methods  for  Three-Dimensional 
Fluid-Structure  Interaction  Problems,"  J.  Acoust  Soc.  Am. ,  Vol .  49,  No.  5, 

(Part  2),  May  1971,  pp.  1505-1510. 

nGeers,  T.  L.,  "Doubly  Asymptotic  Approximations  for  Transient  Motions  of 
Submerqed  Structures,"  0.  Acoust.  Soc.  Am.,  Vol.  64,  No.  5,  Nov  1978,  po.  1500- 
1508. 

12DeRuntz,  0.  A.  and  Geers,  T.  L.,  "Added  Mass  Computation  by  the  Boundary 
Integral  Method,"  Int.  J.  Num.  Meth.  Eng,,  Vol.  12,  1978,  pp.  531-550. 
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FLUID  STRUCTURE  INTERACTION  EQUATION.  Substituting  equation  (2)  into 
equation  (1)  and  equation  (4)  into  equation  (3)  yields  the  coupled  fluid  structure 
interaction  equations. 

Msx  +  C$x  +  K$x  =  fD  -G  Af  (Pj  +  Ps) 

(5) 

Mf  P$  +  pc  Af  Ps  =  pc  Mf  (GTx  -  Oj) 

The  above  equation  (5)  may  be  solved  simultaneously  at  each  time  step  by  the 
transfer  of  -G  Af  P$  and  pc  GTx  to  the  left  side  of  their  respective  equation. 
Such  a  procedure  is  exceedingly  difficult  for  large  systems  because  of  the  large 
coupling  of  the  coefficient  matrices.  Therefore,  a  staggered  solution  procedure 
has  been  developed  that  is  unconditionally  stable  with  respect  to  the  time  step 
for  the  linear  problem13. 

The  computational  strategy  for  the  staggered  solution  procedure  is  embodied 
in  the  following  steps  assuming  the  solution  is  known  at  time  t. 

(1)  Estimate  the  unknown  structural  restoring  force  vector  at  t  +  A  t  from 
the  extrapolation  of  current  and  past  values. 

(2)  Transform  this  extrapolated  value  into  fluid  node  values  and  form  the 
right-hand  side  of  the  fluid  equation,  which  also  involves  the  unknown  incident 
pressure  at  t  +  A  t. 

(3)  Transform  fluid  pressures  into  structural  nodal  forces. 

(4)  Solve  the  structural  equation  for  the  displacement  and  velocity  vectors 
at  t  +  a  t. 

(5)  Transform  the  computed  structural  restoring  force  vector  at  t  +  a  t 
into  fluid  node  values  and  reform  the  right-hand  side  of  the  fluid  equation. 

(6)  Resolve  the  fluid  equation  and  obtain  refined  values  for  the  total 
pressures  at  t  +  a  t. 

(7)  Save  system  response. 

Steps  (1),  (3)  and  (5)  constitute  the  basic  staggered  solution  technique,  while 


13Park,  K.  C.,  Felippa,  C.  A.,  and  DeRuntz,  J.  A.,  "Stabilization  of  Staggered 
Solution  Procedures  for  Fluid-Structure  Interaction  Analysis,"  in  Computational 
Methods  for  Fluid-Structure  Interaction  Problems.  AMD,  Vol .  26,  ASJT7 'Jew  York, 
1?77,  no.  95-124 . 
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Steps  (2)  and  (4)  are  required  because  of  the  difference  between  the  fluid  and 
structural  surface  meshes.  The  iteration  on  the  fluid  solution  reflected  in 
Steps  (6)  and  (7)  has  been  added  to  enhance  accuracy.  Inasmuch  as  the  computation 
time  is  overwhelmed  by  the  structural  solution  requirements,  this  requires  only  a 
small  increase  in  total  run  time.  The  use  of  a  three-step  extrapolation  method 
in  Step  (1)  also  improves  accuracy,  as  discussed  by  Park,  et  al13. 

FLUID  FINITE  ELEMENT  FORMULATION.  By  appropriately  modifying  their  material 
property  matrices,  many  solid  elements,  based  on  the  displacement  formulation,  can 
be  "mocked"  as  fluid  elements.  The  details  of  implementing  this  modification  to 
the  material  matrices  are  given  by  Kalinowski 14.  A  summary  of  the  procedure  used 
is  now  given. 

The  general  relationship  between  stress  and  strain  in  matrix  form  is  given  by 

{c}  =  [E]  {£}  (6) 

where  {0}  and  { e }  are  the  stress  and  strain  vectors  respectively.  The  elasticity 
matrix,  [E],  is  a  6  x  6  isotropic  matrix  which  has  the  general  form 


where  A  and  y  are  the  Lame  constants. 


The  development  of  a  three  degrees  of  freedom  per  node  (3  displacement 
components)  solid  fluid  element  is  straightforward.  By  letting  the  shear  modulus 
(or  equivalently  the  Lame  constant  y)  go  to  zero,  the  dynamic  elasticity  field 
equations15  will  reduce  to  the  displacement  form  of  the  wave  equation 

c2v2u.j  =  u^ ;  c2  =  x/p  (8) 

The  pressure  is  related  to  the  displacement  vector  by 

p  =  -AV-ui  (9) 


14Kalinowski ,  A.  J.,  "Fluid  Structure  Interaction,"  Shock  and  Vibration  Computer 
Proqrams,  edited  by  W.  Pi  1  key  and  B.  Pi  1  key,  The  Shock  and  Vibration  Information 
Center,'  'SVM-10,  1975.  00.  405-452. 

15Kalinowski ,  A.  J.,  "Proceedings  of  the  Fifth  Navy-Nastran  Colloquium,"  Navy 
Nastran  Systems  Office  Report  No.  CMD-32-74,  Seo  1974. 
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where  x  is  the  Lame  constant,  is  the  displacement  vector,  c  is  the  fluid  wave 
speed,  and  p  is  the  fluid  density.  Mote  that  operating  on  equation  (8)  with  the 
divergence  v*()  operator  and  employing  equation  (9)  reduces  equation  (8)  to  the 
pressure  wave  equation 

c2v2p  =  p  (10) 

This  is  the  acoustical  representation  for  the  fluid  field  employed  by  practically 
all  of  the  computer  codes  used  in  fluid  structure  interaction  problems.  The  wave 
equation  (8)  will  be  enforced  by  letting  u  s  0  and  X  =  pc2  in  the  general  stress- 
strain  law,  equation  (7). 

Therefore,  the  local  element  stiffness  and  mass  matrices  for  each  fluid 
element  can  be  created  from  ordinary  solid  linear  elastic  elements  by  appropriately 
choosing  the  values  of  u  and  X.  Elements  created  in  this  manner  have  been  referred 
to  as  "mock"  fluid  finite  elements15.  The  stiffness  matrices  of  the  remaining 
structural  elements  are  then  created  in  the  usual  way.  Finally,  the  total  stiff¬ 
ness  and  mass  matrices  of  the  entire  system  of  elements  are  assembled  in  the  usual 
manner  dictated  by  the  finite  element  method. 

A  boundary  condition  is  also  required  with  fluid  structure  interaction 
methods.  That  is,  at  the  fluid  structure  interface  the  normal  motion  of  the 
structure  equals  the  normal  motion  of  the  fluid.  This  effect  is  achieved  with  the 
mock  elements  by  using  a  double  set  of  nodes  at  the  fluid  structure  interface.  One 
set  belongs  to  the  fluid  and  the  other  to  the  solid.  The  motion  of  the  two  sets 
of  nodes  is  then  assumed  to  be  independent  of  each  other  in  the  tangential 
direction  but  are  forced  to  be  equal  in  the  normal  direction. 

Since  STAGS  does  not  have  a  solid  element,  the  eight  node,  first  order 
hexahedron  element  in  MONSAP16  has  been  modified  according  to  the  procedure 
described  above.  The  mass  and  stiffness  matrices  for  the  mock  fluid  elements, 
generated  using  NONSAP,  are  then  introduced  into  STAGS.  The  USA-STAGS  code  may 
then  be  used  in  the  normal  manner  to  compute  the  response  of  the  fluid-structure 
system. 

USA-STAGS  ANALYSIS 


The  plane  strain  problem  of  the  transient  response  of  two  fluid  coupled 
cylindrical  elastic  shells  to  a  transversely  incident  pressure  pulse  will  serve  as 
a  benchmark  problem  for  the  fluid  finite  element.  The  finite  element  solution  will 
be  compared  with  the  analytical  solution  which  has  been  obtained  by  Huang  for  this 
problem9. 


16Bathe,  K.  J.,  Wilson,  E.  L.,  Iding,  R.  H.,  "NONSAP:  A  Structural  Analysis 
Program  for  Static  and  Dynamic  Response  of  Nonlinear  System,"  Report 
UC  SESM  74-3,  Structural  Engineering  Laboratory,  University  of  California, 
Berkeley,  1974. 
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The  submerged  fluid  coupled  shell  system  and  the  incident  Diane  oressure  wave 
are  depicted  in  Figure  1.  The  fluid  surrounding  the  outer  shell  and  that  between 
the  two  shells  is  considered  to  be  an  ideal  compressible  fluid  and  is  character¬ 
ized  by  its  mass  density,  p,  and  sound  speed,  c.  The  Doubly  Asymptotic  Approxima¬ 
tion  (DAA)  is  used  for  analyzing  the  transient  interaction  of  the  infinite  fluid 
exterior  to  the  outer  shell,  the  incident  pressure  wave,  and  the  outer  shell 
itself.  This  leaves  only  the  inner  and  cuter  shells  and  the  entrained  fluid  tc  be 
modeled  using  finite  elements.  Due  to  symmetry,  modeling  only  the  region  from 
0<6<tt  i$  necessary.  The  finite  element  model  used  is  shown  in  Figure  2.  It 
contains  36  plate  bending  elements  for  the  inner  and  outer  shell  and  72  hexahearcr 
fluid  elements  for  the  entrained  fluid. 

The  middle  surface  radii,  thickness,  mass  densities.  Young’s  moduli  and 
Poisson's  ratios  of  the  outer  and  inner  shells  are  (ae,  he,  p®,  Ee,  \>e)  and 
(a,  h,  p$,  E,  v),  respectively.  The  material  properties  and  dimensions  used  in 
the  checkout  problem  are: 

p  =  pe  =  0.C362  lb/in3 
c  =  ce  =  60C00  in/sec 
Pes  =  Ps  =  0.282  lb/in3 
Ee  =  E  =  30  x  106  psi 
\>e  s  v  =0.3 
ae  *  6.25  in 
a  =  5.0  in 
he  =  0.0363125  in 
h  =  0.14525  in 

RESULTS 

An  analysis  of  the  submerged  fluid-coupled  shell  was  done  for  the  case  cf  a 
step  wave  of  infinite  duration.  The  USA-STAGS  calculations  were  carried  out  for 
600  steps  using  a  time  step  of  about  1/2C0  of  the  transit  time  (the  time  for  the 
incident  wave  front  to  traverse  the  distance  of  one  outer  shell  diameter). 

Figure  3  shows  the  time  history  of  the  radial  velocity  (w)  at  various  locations 
of  the  inner  shell.  (Note  that  the  velocity  has  been  non-dim, ensicr.alized  by 
ft=P0/pc,  where  pQ  is  the  magnitude  of  the  incident  pressure  wave.  The  non- 

dimensional  time,  T1 ,  is  given  by  T'*(T-n)/c  where  T=cet/ae,  na(l-s),  q=a/ae,  and 

t  is  the  actual  problem  time.  Physically,  the  motion  of  the  outer  shell  starts 

at  T=0,  while  that  of  the  inner  shell  starts  at  T=n=(l-a/ae) . )  In  general,  the 
agreement  between  the  USA-STAGS  solution  and  Huang's  analytical  solution  is  quite 
good  though  the  calculated  results  contain  seme  higher  frequency  oscillations. 

The  late  time  translational  velocities  of  the  shell  compare  very  well. 

A  comparison  of  the  analytical  and  USA-STAGS  solutions  for  the  time  histories 
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of  the  hoop  stress  resultants  at  various  locations  of  the  Inner  shell  Is  shown  In 
Figure  4.  (Note  that  the  stress  resultants,  N0,  have  been  non-dimenslonallzed  by 

fl=P0/pc2.  The  time,  T' ,  is  non-dimenslonallzed  as  before.)  Again,  the  agreement 

between  the  two  solutions  is  quite  good.  Also  shown  on  this  figure  are  results 
from  Huang's  analytical  solution  done  without  the  presence  of  an  outer  shell. 

CONCLUSION 


Comparison  of  the  analytical  and  USA-STAGS  solutions  to  the  problem  of  two 
submerged  fluid-coupled  shells  Impinged  upon  by  a  transversely  incident  plane 
pressure  wave  is  quite  favorable.  It  is  apparent  from  this  analysis  that  the  use 
of  the  "mock"  fluid  finite  element  is  an  effective  technique  for  solving  this  class 
of  fluid-structure  interaction  problems. 


/ 
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FIGURE  2  FINITE  ELEMENT  MODEL  FOR  FLUID-COUPLED  SHELL  SYSTEM 
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FIGURE  3  TIME  HISTORIES  OF  RADIAL  VELOCITIES  OF  THE  INNER  SHELL 
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FIGURE  4  TIME  HISTORIES  OF  THE  HOOP  STRESSES  AT  THE  MIDDLE  SURFACE 
OF  THE  INNER  SHELL 
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